/*
 * Matrix.h
 *
 *  Created on: Apr 9, 2010
 *      Author: andrew
 */

#ifndef MATRIX_H_
#define MATRIX_H_


#include "Vector1D.h"
/**
 * uses fortran BLAS and LAPACK libraries
 */
extern "C"{

// matrix vector operations
void dgemv_(const char *trans, const int *m, const int *n, const double *alpha,
		const double *A, const int *lda, const double *X, const int *incx, const double *beta,
		double *Y, const int *incy);

// solve system of linear equations
void dgesv_(const int *n, const int *nrhs, double *A, const int *lda, int *ipiv,
		double *B, const int *ldb, const int *info);

// matrix matrix multiply
void dgemm_(const char *transa, const char *transb, const int *m, const int *n,
		const int *k, const double *alpha, const double *A, const int *lda, const double *B,
		const int *ldb, const double *beta, double *C, const int *ldc);

// matrix exponential
void dgpadm_(const int *ideg, const int *m, const double *t, const double *H,
		const int *ldh, double *wsp, const int *lwsp, const int *ipiv, int *iexph,
		const int *ns, const int *iflag);

// cholesky factorization
void dpotrf_(const char *uplo, const int *n, double *A, const int *lda, const int *info);

// LU decomoposition of a general matrix
void dgetrf_(int* M, int *N, double* A, int* lda, int* ipiv, int* info);

// generate inverse of a matrix given its LU decomposition
void dgetri_(int* N, double* A, int* lda, int* IPIV, double* WORK, int* lwork, int* info);

void dger_(const int *m, const int *n, const double *alpha ,double *X,
		const int *incx,double *Y, const int *incy, double *A, const int *lda);

void dgels_(const char *transa, const int *m, const int *n, const int *nrhs,
		double* A, const int *lda, double *b, const int *ldb, double *work,
		const int *lwork, const int *info);

void dgtsv_( const int *n, const int *nhrs, double *dl, double *d, double *du,
		double *b, const int *ldb, const int *info );
}

class Matrix {
	double* matrix;
	int m,n;

public:
	Matrix();
	Matrix(int rows, int cols);
	virtual ~Matrix();
	Matrix(const Matrix& copyme);
	Matrix& operator=(const Matrix &mObj);

	double& operator()(int i, int j) const;
	Vector1D operator*(const Vector1D& right) const;
	Matrix operator*(const Matrix& B) const;
	Matrix operator*(const double s) const;
	Matrix operator+(const Matrix& B) const;
	Matrix operator-(const Matrix& B) const;
	Vector1D solve(const Vector1D& b) const;
	Matrix solve(const Matrix& B) const;
	Matrix getTranspose() const;
	Matrix expm(double tstep) const;
	Matrix appendMatrixBelow(Matrix& A) const;
	Matrix appendMatrixSide(Matrix& A) const;
	Matrix appendVectorCol(Vector1D& v, int col) const;
	Vector1D selectCol(int col) const;
	Matrix initConst(double init) const;
	Matrix diag(Vector1D d) const;
	Matrix subselect(int rowStart, int colStart, int nrows, int ncols) const;
	Matrix blkdiag(const Matrix& B) const;
	Matrix chol() const;
	Matrix inv() const;
	static Matrix eye(int n);
	double rows(){ return m; }
	double cols(){ return n; }
	Matrix selectEntries(const Vector1D& selectRows, const Vector1D &selectCols);
	void printMatrix(void) const;
	void printMatrix(int rowstart, int rowend, int colstart, int colend) const;
	Vector1D vecNot(int i, int j, int k) const;

	static Matrix fortran2c(const double *B, int n, int m);

private:

	static double* c2fortran(const Matrix A);

};

#endif /* MATRIX_H_ */
